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Abstract. The formation and evolution of galactic disks 
are complex phenomena, where gas and star dynamics are 
coupled through star formation and the related feedback. 
The physical processes are so numerous and intricate that 
numerical models focus, in general, on one or a few of them 
only. We propose here a numerical model with particu- 
lar attention to the multiphase nature of the interstellar 
medium; we consider a warm gas phase (> lO'^K), treated 
as a continuous fluid by an SPH algorithm, and a cold gas 
phase (down to lOK), fragmented in clouds, treated by a 
low-dissipation sticky particles component. The two gas 
phases do not have the same dynamics, nor the same spa- 
tial distribution. In addition to gravity, they are coupled 
through mass exchanges due to heating/coohng processes, 
and supernovae feedback. Stars form out of the cold phase, 
and re-inject mass to the warm phase through SN explo- 
sions and stellar winds. The baryons are embedded in a 
live cold dark matter component. 

Baryonic disks, initially composed of pure gas, encounter 
violent instabilities, and a rapid phase of star formation, 
that slows down exponentially. Stars form in big clumps, 
that accumulate in the center to build a bulge. Exponen- 
tial metallicity gradients are obtained. External infall of 
gas should be included to maintain a star formation rate 
in the disk comparable to what is observed in present disk 
galaxies. 



1. Introduction 

The formation of a galactic disk involves such complex 
physical processes, that only numerical computations can 
give some insight in their relative importance. The essen- 
tial force is gravity, which is nowadays easily controlled 
through N-body codes, based either on grid calculations 
(such as EFT), or on hierarchical grouping of particles 
(such as the tree-code). Both algorithms optimize the com- 
puting time, growing roughly as A'^ln(iV). The main im- 
provement in this domain is obtained by going towards 
higher and higher spatial resolution, with an ever growing 
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number of particles. The gravitational force is softened at 
small scale, to reduce unphysical two-body relaxation, re- 
sulting in a resolution equal or slightly smaller than the 



inter-particle distance (Romeo 199^, Dehnen 2001, Knebe 



at al. 2001) 



It has been recognized for a long time that the dis- 
sipative component is also an essential feature in galaxy 
evolution, and the interstellar medium has been widely 
introduced in galaxy simulations, either as a continu- 
ous fluid (Eulerian grid codes, van Albada & Roberts 



1981, Lagrangian SPH codes, Hernquist and Katz 1989) 



or through sticky particles algorithms (Combes & Gerin 



1985). However, contrary to the pure gravity case, the re- 



sult of simulations now depends not only on the spatial 
resolution and computer power, but on the physical as- 
sumptions made on the nature of the ISM and related 
processes. 

Star formation and feedback, coupling the stars and 
gas by non-gravitational processes are essential (Katz 



1992, Mihos & Hernquist 1994). Since the detailed pro- 



cesses involved in star formation are not yet well known, 
this introduces further uncertainties and liberties in the 
modeling. The most widely adopted recipe to control star 
formation is based on a local Schmidt law (i.e. the star 
formation rate is locally proportional to some power of 
the volume density), sometimes associated with a thresh- 
old for star formation (e.g. Friedli & Benz 1995). However, 
the local Schmidt law is not actually observed in galaxies, 
and the justification comes from an observed global em- 
pirical Schmidt law, averaged over the whole galaxy (e.g. 



Kennicutt 1998). Other star formation recipe s are based 



on Jeans inst abilit y (e.g. Steinmetz & Miiller 1994 Ger- 
ritsen & Icke 1997 ), or cloud-cloud collisions (Noguchi & 
Ishibashi |l986| ). 

Due to the limited number of particles, each star "par- 
ticle" represents in fact a stellar cluster of the order of 
lO'^ to 10^ M0. The conversion of a fraction of a gas par- 
ticle into stars requires the creation of a large number of 



star-particles (e.g. Katz 1992), or the consideration of "hy- 
brid" particles, that are transiently containing both gas 
and stars, and have the same dynamics for a while (Mihos 



Hernquist 1994). Alternatively, some "starlet" parti- 
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cles are transiently created, and merged with the nearest 
neighbors (Jungwiert et al. 2001). Whatever the choice, 
some approximations are made, related to the limited res- 
olution in time, mass and space, respectively of the order 
of 10^ yr, 10^ M0 and 0.3 kpc for a typical giant spiral 
galaxy. 

The same limitations occur when mass loss and feed- 
back are considered. Eith er on ly stellar heating is con- 
sidered (Gerritsen & Icke 1997), or the supernovae me- 



chanical energy only(Mihos & Hernquist 1994), or more 
sophisticated models are used including several types of 
supernovae, planetary nebulae, stellar winds, evaporation 
and condensation (e.g. Theis et al. 1992, Samland et al. 



1997, Berczik 1999). 



Stellar nucleosynthesis can be followed to describe the 
detail ed me tal enrichment history of the galaxy (e.g. Lia 
et al. 2001). Models reveal that the frequently used in- 
stantaneous recycling is a too simple approximation (see 
also Jungwiert et al. |200l| ). 

Of first importance is the thermal evolution of the gas, 
since thermal instabilities are responsible for cloud con- 
densation, feeding the cold ga s pha se, and thus star for- 
mation (Hultman & Pharasyn 1999| ). Cooling and heating 
processes are taken into account depending on the tem- 
perature, density and metallicity of the g as (e. g. Dalgarno 
& McCray |l972| , Sutherland & Dopita [19931) . However, 
the time-scales involved in the thermal processes can be 
much smaller than the dynamical time-scales, and the spa- 
tial scales of the corresponding processes are much below 
the spatial resolution of the simulations. Phenomenolog- 
ical recipes are therefore used to convey at large-scales 
the resulting effects of the processes occurring below the 



resolution (e.g. Thacker & Couchman 2000): the heating 
energy is smoothed over the resolution scale, and the ef- 
fective time-scales slowed down. 

Many simulation models consider only one gas phase, 
treated as a continuous fluid at the virial temperature cor- 
responding to a galaxy potential (~ 10^ K). A simple ap- 
proximation is a strictly isothermal gas, since the cooling 
is very efficient above 10"* K (Barnes & Hernquist 1991, 
Mihos & Hernquist 1994). Some allow gas to cool down 
to low temperatures (~ 10 K), and to spread into several 



temperature phases (e.g. Gerritsen & Icke 1997; Yepes et 
al. 1997, Thacker & Couchman 2001). However, their dy- 
namics is still that of a one-phase medium, given the in- 
sufficient density contrast that can be dealt with the nu- 
merical simulations on large-scales. Significant contrasts 
can be achieved only in simulations of small volumes, that 
approach more realistically the multiphase nature of the 
interstellar medium (e.g . Rosen & Bregman ( 1995 ), Wada 
& Norman ( |l999i ^00l|) . 

In the present paper, we implement most of the above 
processes to investigate the formation and evolution of 
disk galaxies (gravity, stellar and gas dynamics, star for- 
mation, feedback and massloss heating and cooling, metal 
enrichment..), and we focus on the large-scale multi-phase 



dynamics of the gas. The "warm gas" component is mod- 
eled by a continuous fiuid with an SPH code, and contains 
a large range of temperatures, up to hot gas re-injected by 
supernovae, though the bulk of the gas is found between 
lO'' and 2 lO'^K. The "cold gas" component corresponds to 
the cloudy and fragmented, essentially molecular medium, 
between lOK and lOOK. It is modeled as a separate phase, 
via individual cloud-particles, subject to inelastic colli- 
sions. Since at large galactic scales, the density contrast 
necessary to form such a fragmented structure (12 or- 
ders of magnitude) is far beyond the present computa- 
tional power, the cloud-particle approximation appears 
well suited. For the cloudy component, the pressure forces 
or viscosity forces are negligible. The dissipation occur- 
ring during cloud collisions at AU scales, much below our 
spatial resolution, is modeled phenomenologically. 

Two general classes of models have been proposed for 
the formation of galaxies: either big galaxies can form 
through the monolithic collapse of a very massive bary- 
onic system (e.g. Eggen, Lynden-Bell & Sandage 1962), or 
dwarf galaxies form first, and act as the building blocks of 
larger systems, that subsequently form through recursive 
mergers, in the hierarchical scenario (e.g. Kauffmann et 
al 1999). The first scenario assumes that the gas experi- 
ences violent 3D star formation, so quickly that it has no 
time to settle into a disk. This has been proposed to form 
large elliptical galaxies or bulges. However, it is today rec- 
ognized that many (if not all) elliptical galaxies can form 
through mergers of spiral galaxies, or smaller entities, as 
first suggested by Toomre & Toomre (1972), and devel- 
oped further by Schweizer (1986), while bulges can form 
similarly by mergers, and also through secular evolution 
(e.g. Combes 2000). Although the most successful cosmo- 
logical scheme today, based on inflation, promotes primor- 
dial density fluctuations that are scale-invariant, which to- 
gether with the existence of cold dark matter, results in hi- 
erarchical formation of structures, the collapse of baryons 
inside dark matter halos is still an unsolved problem, with 
many uncertainties and free parameters. The rate at which 
stars are formed in a given system, and whether the stellar 
populations observed today can be attributed to a sudden 
event like the monolithic collapse, or a more progressive 
one like the gradual and recursive merging, is a matter 
of debate. It is likely that both point of views may corre- 
spond to some observed galaxy examples today, since the 
major merger of two gas rich giant spirals have some sim- 
ilarities with a monolithic event. More subtle diagnoses 
should be searched for, and large number statistics should 
be used to determine the dominant processes. In this work, 
we focus on the formation and evolution of a single galaxy. 
To achieved the required resolution we use the monolithic 
collapse scenario. 

In section 2 we detailed the assumptions of the multi- 
phase model, and present the code and numerical methods 
in section 3. In section 4, 5 and 6 we present and discuss 
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Fig. 1. Schematic presentation of the model. The four 
phases are represented with the main physical processes 
responsible for mass and energy exchanges. 

the results from our simulations. We give our conclusions 
in section 7. 



2. The multiphase model 

We will describe in detail the different aspects of the 
model. Fig. 1 presents a schematic view for quick refer- 
ence. 



2.1. Description of the four components 

2.1.1. Dark matter 

We have compared the results of two usual methods to 
take dark matter into account. The first uses an analytic, 
static potential. The second uses collisionless particles. 
The first choice obviously has a major shortcoming: the 
dark matter density field does not respond to the evolution 
of the baryonic matter. The static density configuration of 
the dark matter affects the evolution of the system, and 
phenomena such as transfer of angular momentum from 
baryonic matter to dark matter are ignored. While using a 
static potential may be considered in the case we focus on 
in this work; the evolution of a single galaxy within a dark 
matter halo, it would appear unreasonable in the case of 
mergers where dark matter halos undergo large deforma- 
tions. We will show in this paper that some difficulties also 
appear when using collisionless dark matter particles with 
insufficient numerical resolution. We show in Sec. 5 that 
using a "small" number of heavy dark matter particles 
results in unphysical heating of the stellar disk within 2 
Gyr. This problem can be overcome partly by using large 
smoothing lengths for the dark matter, but more reliably 
by increasing the number of particles. 

2.1.2. Stars 

Stars form the second type of collisionless matter. They 
interact with the gas components through processes de- 



scribed in Sec. 2.2 . In this work, the mass of each stellar 
particle is of the order of ^ IO^Mq in most simulations 
and ~ 10^ Mq in the high resolution simulation. Each par- 
ticle thus represents a large number of stars, and such 
properties as supernovae rate or metallic enrichment will 
depend on the assumed initial mass function (see Lia et 
al. 2001). In this work, we do not study the detailed de- 



pendence on the IMF, but rather use standard values for 
quantities such as supernovae rate. Our model includes 
many processes and we choose a simple representation for 
each of them. 



2.1.3. Warm gas 

The warm gas in the ISM is a highly complex system. 
Dynamically, it can be described in first approximation as 
a self-gravitating dissipative fluid. In our model, the warm 
gas phase is composed of SPH particles obeying the ideal 
gas equation of state: 

P^{^~l)pu (7 = ^). (1) 

Here, P is the pressure, p the density, and u the thermal 
energy per unit mass. The thermal energy is governed by 
equation: 



du ~ 2 A / ^ 

p— + PVv = TiffFuv + nfjA{u), 



(2) 

where v is the velocity of the gas and uh the number den- 
sity of hydrogen. P is the viscosity-corrected pressure (see 
Sec. 3.1 for details and references on the SPH implemen- 
tation). On the right hand side are heating and cooling 
terms described in detail in Sec. 2.2.1. The evolution of 
the internal energy is limited to the corresponding finite 
range in temperature llOOOK -10^ K. The lower limit cor- 
responds to the transition to the cold star-forming gas 
phase. The upper limit is for numerical convenience, since 
in this work most SPH particles lie in the 11000 K - 20000 
K range, at the virial temperature of typical galaxies. 

2.1.4. Cold gas 

A number of numerical models stop the cooling of SPH 
particle at 10000 K (e.g., Weil et al |1994 Thacker & 



Couchman 200C). This can only be a first approximation. 

Indeed, a large amount of the gas in galaxies is ob- 
served in cold ( 10-100 K ) dense molecular clouds. These 
objects are essential to the global dynamics of the galaxies 
since they are the location of star formation. Their phys- 
ical properties and dynamics differ strongly from those of 
the surrounding diffuse warm gas. Indeed, each molecular 
cloud is in local virial equilibrium and the overall filling 
factor of the cold gas is very small. Consequently, it can- 
not be accurately described as a fiuid (which supposes a 
rather continuous medium). Moreover, current treatment 
of SPH uses the ideal gas state equation to relate pressure, 
density and internal energy of particles which, at present 
days numerical resolution, have the mass of a giant molec- 
ular cloud. This ideal gas description is hardly adequate 



4 



Semelin & Combes: Formation and evolution of galactic disks with a multiphase numerical model 



for cold, self-gravitating, GMC-like particles, which should 
rather obey Larson laws (Larson 1981). For these reasons, 
we feel that a specific non-SPH treatment should be given 
to these star-forming gas particles. 

Some attempts have been made in this direction. Let 
us first mention the works by Wada & Norman (1999) and 
Wada (2001) who include proper cooling down to 10 K in 
Eulerian simulations of small regions of the galaxy ( ~ 100 
pc). Yepes et al. ( 1997[) , using a PPM scheme, and Hult- 
man & Pharasyn ll999| ), using a SPH simulation, both 
identify a phase of cold, star-forming gas. However, both 
still treat the dynamics of the gas as a single fluid. Relying 
on the fact that the largest fraction of the gas in the ISM is 
in the form of cold clumpy fragments (molecular clouds), 
Noguchi (1999) adopts a sticky particle scheme to describe 
the coUisional gas. He does not account for the warm dif- 



fuse gas. Andersen & Burkert (2000) develop a more elab- 



orate model to treat the cold gas clouds as a coUisional 
fluid. They take warm gas into account only through a 
constant external pressure exerted on cold clouds. In this 
work we attempt to fully model a cold and a warm gas 
phase with different dynamics. 

While we describe the warm gas with the usual SPH 
technique, the cold gas particles experience gravitational 
forces only, and are subject to inelastic collisions follow- 



ing the sticky particle scheme (Levinson & Roberts 1981). 
Within one time step, collisions occur with a probability 
P ~ Aoc^v, where pi^^ is the local cold gas density, and 6v 
is the local velocity dispersion. The only "thermal" evo- 
lution for cold gas particles is the possibility to be evap- 
orated back to the warm gas phase through supernovae 
feedback (see Sec. 2.2.3) 

2.2. Physical processes 
2.2.1. Cooling and heating 

The evolution of the thermal energy of warm gas particles 
follows equation (2) . The sources of variation of the ther- 
mal energy are the following: adiabatic expansion, heat- 
ing from viscosity, heating from an external radiation field 
(the Fuv term) and radiative cooling (the A term). 

We use the tabulated values of A in Sutherland & Do- 
pita ( 1993|) . We have not included the effect of metaUic- 
ity on the cooling: we use a constant solar metallicity in 
the cooling function. Indeed, most of the gas lies in the 
10000K-20000K range, where metallicity has little effect. 
In the typical conditions of galactic dynamics, the cooling 
time can become much shorter than the dynamical time in 
resolved dense regions. Reducing the dynamical time step 
accordingly is not appealing since it would drastically in- 
crease de CPU requirements. Moreover, it would not be 
consistent with the time scale resolution of gravitational 
effects. Consequently we choose to damp the thermal fluc- 
tuations (the same solution is applied by Weil et al. 1998] ) . 
Whenever the variation of the internal energy of a particle 



within one time step is larger than one fourth of its current 
value, it is damped to one fourth. In addition, adaptative 
small steps are used to follow the very fast variations of 
A as a function of u in the 10000K-20000K range (this is 
similar to the integral scheme use by Thomas & Couch- 
man 1992). This results in a stable thermal behavior for 



-24 



the gas. 

A constant uniform heating F^v ^ 10~" erg s^" is 
applied to the gas to model the background UV radiation 
field. The variations of F only shifts the mass balance be- 
tween warm and cold phase as other parameter can, such 
as the minimal temperature of the warm gas. It is how- 
ever important to include a non-zero F to help sustain the 
warm phase. In our study we keep F constant. Gerritsen 



& Icke (1997) propose a detailed model for the heating. 



The second source of heating comes from supernovae 
feedback. This is the only way in our model for cold gas to 
return directly to the warm gas state. It will be described 
in Sec. 2.2.3. 

Let us explained now how the transition from warm 
gas to cold gas is handled. A parameter uq is set which 
defines the thermal energy at the transition. It is a free 
parameter, which value corresponds to a temperature of 
11000 K. A simple method would be to transform into cold 
gas any warm gas particle cooling down to uq. However, 
owing to the different density dependence of cooling and 
heating terms in eq (2), this choice induces a sharp density 
threshold for turning warm into cold gas. When the aver- 
age warm gas density falls noticeably below the threshold 
the transfer to cold gas is stopped and the cold gas phase 
is progressively depleted by star formation. We obtained 
better equilibrium between the phases with a smoother 
transition. Accordingly, and for the purpose of transition 
to cold gas only, the thermal energy of each SPH parti- 
cle of warm gas is considered as a Gaussian distribution 
around a central value u with a dispersion of au, rather 
than a Dirac distribution d{u). This makes sense from a 
physical point of view since the scale resolution of the sim- 
ulation is about 100 pc, leaving out thermal fluctuations 
below this scale, which can be crudely taken into account 
using the Gaussian distribution. We use a cr„ correspond- 
ing to ctt ^ 1500K. The tail of the Gaussian distribution 
extending below uq can be considered as cold gas. Conse- 
quently we can define the fraction fc of cold gas in a warm 
particle as 



fc = 0.5(1 - erf(- 



Uq 



))■ 



This is the fractional mass that is used to produce cold 
gas particles following the algorithm described in Sec. 3.2. 

2.2.2. Star formation 

There are two steps involved in modeling star formation. 
First, select under which conditions stars can be formed 
from the gas, then specify at which rate the formation 
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occurs. Katz (1992) proposes a set of criteria; stars can 
form only in Jeans unstable regions and within a conver- 
gent flow of the gas. In these regions they forms at a rate 
following a Schmidt law. Navarro & White ( 1993| ) add the 
condition that the star forming regions are cooling rapidly. 
They reduce the cooling time criterion and the Jeans in- 
stability criterion to a density threshold. This is coherent 
with results derived from the observations (Martin & Ken- 
nicutt 2001). Other authors argue that the limited reso- 
lution of the simulations, where a gas particle has usually 
the mass of a whole GMC, makes the use of dynamical cri- 
teria unreliable (e.g Mihos & Hernquist 1994, Hultman & 



2.2.4. Stellar mass-loss 



Pharasyn 1999). We agree with these considerations and 



use a simple Schmidt law to describe star formation: 



dt 



(3) 



C is a constant combining a star formation efficiency, a 
typical formation time and a typically density. All cold 
gas particles are candidate for star formation and the star 
formation rate is fixed by a Schmidt law with index 1.5 
or 1. Finally let us mention the work by Gerritsen & Icke 
(Gerritsen & Icke 1997 ), who solely rely on dynamical cri- 
teria such as convergent flow and Jeans instability, and 
attempt to derive a Schmidt law as a result of the simu- 
lations. 



2.2.3. Supernovae feedback 

A lower estimate for the amount of energy released in su- 
pernovae for a Salpeter IMF is 10^^ erg per solar mass of 
formed stars. During the first 2 Gyr of its life, a galaxy 
converts a large fraction of its gas content into stars. The 
total energy released by supernovae during this process 
is many times larger than the total kinetic energy of the 
galaxy. It appears important to take this phenomenon into 
account. However since most of the feedback is in the form 
of thermal energy in the dense gas environment of young 
stars, it will be quickly radiated away. A small fraction 
is also fed back as kinetic energy in expanding gas shells. 
This phenomenon is strong enough to eject gas from dwarf 
galaxies. For a a detailed analys is of the different algo- 
rithms see Thacker & Couchman (2000). In our model, we 
include both thermal and kinetic feedback. Kinetic energy 
is imparted to any neighboring gas particles in the form 
of a velocity kick opposite to the forming star direction. 
Thermal energy is directly deposited on a small number of 
neighboring cold gas particles ('^ 3), which are converted 
into warm gas particle and brought to a temperature of 
20 OOOK to 100 OOOK depending on the simulation. The 
thermal feedback is largely affected by the damping of the 
thermal evolution of SPH particles described in Sec. 2.2.1. 
It actually has an effect similar to imposing an adiabatic 
period of a few Myr as described in Thacker & Couchman 
(2000). We will show in Sec. 6.1 that the thermal feedback 



Jungwiert et al. ( ^001 ) show that within 10 Gyr, super- 
novae explosions and stellar winds can return over 40 % 
of the stellar mass to the ISM for a star population fol- 



lowing Scab's IMF (Scalo |199^ ). Of these 40 %, less than 
10 % are expelled in supernovae explosions. Conseq uently , 
the instantaneous recycling approximation (Tinsley 1980 ) 
used in most numerical simulations including star forma- 
tion and metallicity enrichment (e.g. Hultman & Pharasin 



19991) can only be a first step. Indeed, stars keep releasing 



associated to the damping has a definite effect on the mass 
balance between the phases. 



gas through stellar winds long after they have separated 
from the gas clouds where they were born, and refill the 
gaseous disk proportionally to the smooth star density 
field, not onl y in lo cations where stars are formed. Rosen 
& Bregman (1995) improved on the instantaneous recy- 
cling approximation by using a constant mass loss rate 
for the stars. Jungwiert et al. ( ^001 ) and Lia et al. (2001) 
propose much more detailed models. They derive a time 
dependent mass loss rates for various IMF, and include 
them in simulations. In our model we use a simple version 
of Jungwiert et al. mass loss rate: 

^ ^ Af^^„ ^ . (4) 

dt {t - tbirth + To) 

We use the constant values c — 0.055 and T — 5 Myr. Our 
algorithm for mass exchange makes it unnecessary to deal 
the with merging of small star particles as in Jungwiert et 
al. (see Sec. 3.2). 

2.2.5. Metal enrichment 

The production of heavy elements in stars and subsequent 
ejection in SN explosions or stellar winds lead to the mod- 
ification of the chemical composition and cooling proper- 
ties of the ISM. Since our model includes both stellar for- 
mation and stellar mass loss, it is natural to follow the 
chemical evolution of the ISM. However, as mentioned in 
Sec 2.2.1, we will not take the modifications of the cooling 
properties into account. 

We use the yield y, defined as the mass of metal pro- 
duced per unit mass of stellar hydrogen (Tinsley 198C). If 
a star was initially formed from gas with a metallicity Zs , 
the gas produced by this star though the mass loss process 
has a metallicity: 

Zg = Zs + y{l - Zs) 

During the gathering process used to produce gas par- 
ticles or stars (see Sec. 3.2), the metallicity of the different 
contributions is averaged. This introduces a numerical dif- 
fusion at the resolution scale. The physical processes lead- 
ing to metal enrichment of the gas are stellar winds and 
supernovae explosions. The diffusion of metallicity in the 
actual ISM is probably realized through turbulent convec- 
tive mixing. 

The main difference of our model with the instan- 



taneous recycling approximation (Tinsley 1980) is that 
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it works within the continuous mass-loss scheme, so en- 
riched gas is not only released in star forming regions but 
throughout the stellar disk. 

3. Numerical methods 

3.1. Tree-SPH 

We use a N-body representation of the system. Gravita- 
tion and hydrodynamics are the main factors driving the 
evolution. To compute the corresponding forces acting on 
each particle, we use the now familiar tree-SPH method 
(Hernquist & Katz |l989[ ). 

Gravitation forces are co mputed using the so-called 
tree algorithm (Barnes & Hut 1986 ). This particular im- 
plementation was u sed in previous works, for example 
Semelin & Combes ( 2000 ). Following this algorithm, the 
contribution by a group of distant particles to the grav- 
itational force acting a given particle can be computed 
using the multipole moments of the group mass distribu- 
tion. It is actually computed if the angular size of the 
group, seen from the particle, is smaller than an opening 
criterion 9. Otherwise, the group is divided into subcom- 
ponents. Throughout this work we use a value of 6* = 0.8. 
The multipole expansions are carried out to quadrupolc 
moments. With this setup, the typical error on the gravi- 
tational force is smaller than 1% (Hernquist |1987 ). 

Hydrodynamic forces are computed through the SPH 
method (Gingold & Monaghan |l977| , Lucy |l977| ). For a 
complete description of the method, see, for example, 
Monaghan (1992). However, implementations of the SPH 
technique come in many different flavors (see Thacker et. 
al. (2000) for a comparative study). Let us summarize 
our choices for this implementation. We use the ideal gas 
equation of state: 

P=(7-l)pu (7 = ^) 
We use a kernel based on a spherically symmetric spline 



function (Monaghan & Latanzio 1985). The individual 
smoothing length of each particle is adaptative and we 
use the arithmetic average hij = !h±!h. compute the 
interaction between two particles. We have taken care to 
insure that whenever particle B acts as a SPH-neighbor 
of particle A, then the reverse is true, particle A acts 
as a SPH-neighbor of particle B. This is important for 
a good energy conservation. The time update of the in- 
dividual smoothing le ngth follows the scheme proposed 
by Hernquist & Katz (1989). Finally, we use the viscosity 
described in Monaghan ( 1992 ). 

To measure the quality of our implementation, we run 
a simulation that is becoming a standard validation test 
for this type of codes: the collapse of an initially static, 
isothermal sphere of self-gravitating gas. This was first 
studied by Evrard ( 198^ ) and has been reproduced by 
many authors since then ( e. g., Hernquist & Katz 1989, 
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Fig. 2. Evolution of total, thermal, potential and kinetic 
energies of an initially isothermal sphere of gas with den- 
sity profile p{r) ~ i. We use dimensionless unit, and 
G = 1. The thermal and potential extrema correspond 
to the instant of maximal collapse. 



Thacker et al. 200C, Springel et al. 2001,..). Let us consider 
a sphere of radius i?, mass M, density profile 



M 1 



27ri?2 r ' 



and uniform internal energy 



0.05- 



.GM 



The sphere is composed of TV particles initially at rest, 
distributed within R with the proper probability to realize 
the profile p(r). We use R — 1, Al — 1 and G — 1. Since 
our implementation uses a single time step, we have to 
choose a small value. We set St — 0.001, and integrate 
from i = to t = 3. We use TV = 10000 particles, a 
softening e = 0.02 and an opening criterion 6 = 0.8. The 
evolutions of total, kinetic, potential and internal energies 
are given on Fig 2. The maximum relative variation of 
the total energy is 0.3 %. The evolutions on Fig 2. match 
closely those given by Thacker et al. (2000) or Springel et 
al. (2001). This result was obtained using 40 neighbors for 
each SPH particle. We have observed that using a smaller 
number of neighbors, 15, for a constant TV, produces a 
similar evolution with a thermal peak 35 % higher, just 
as larger N produce higher peaks. This boils down to the 
degree of smoothing of the central collapse by the SPH 
algorithm, the smaller the smoothing length, the larger 
the thermal peak. 
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3. 2. Implementation of mass exchanges 

As presented in sec. 2, our model takes into account four 
different types of particles: dark matter, stars, warm gas 
and cold gas. As a result from cooling, warm gas can be 
converted to cold gas. In dense regions, cold gas is then 
converted to star following the star formation rate fixed 
by a Schmidt law. In turn, stars expel warm gas following 
the mass- loss prescription (eq. 4). And finally, cold gas can 
be transformed into warm gas as a result of supernovae 
feedback. Except for the vaporization of cold gas to warm 
gas, which receive a simple treatment, we use a common 
numerical procedure to realize these exchanges. 

A number of procedures have been used in the liter- 
ature, applied mostly to star formation. In one way or 
another, all procedures use a variable to keep track of 
the fractional mass of stars already produced in a given 
gas particle (this can be generalized to any two phases). 
We can then separate procedur es in two types. In the 
first type, pioneered by Katz ( 1992| ), the gas particle 
spawns small star particles each time its fractional star- 
mass reaches a chosen low threshold value. In this proce- 
dure the dynamics of coUisionless star particles is differen- 
tiated early-on from the dynamics of the gas. The draw- 
back is that it leads to the creation of many star particles 
with small masses. This increases the CPU cost and leads 
to two-body relaxation problems: the young stellar disk, 
composed of low-mass particles will be sensitive to the dy- 
namical heating by heavier particles. We will show in sec. 
5 that this is a sensitive issue, especially with heavy dark- 
matter particles. Jungwiert et al. (2001) have applied this 
procedure to the mass-loss of the stars. They have taken 
care of some of the drawbacks by including schemes to 
merge small-mass particles of the same type. 

In the second type of procedure, the entire gas par- 
ticle is converted to star when the fractional star-mass 
reaches a high threshold, 80% for example. This was used 
by Mihos & Hernquist ( 1994 ). They also choose to sub- 
tract the fractional star-mass from SPH density evalua- 
tions to avoid having a 70 % star particle behave as a full 
gas particle. Oth er authors choose not to, e.g. Thacker & 
Couchman (2000). An extension of this procedure is to use 
a probabilistic approach to convert the gas particle into a 
star particle, with a probability function depen ding on the 
fractional star-mass (see, e.g., Lia et al. 2001 ). The main 
drawback is that large amounts of star-mass may have to 
follow gas-particle trajectories for a long time. 

Our scheme, similar to Steinmetz & MuUer ( |1994[ ), 
combines the advantages of keeping the number of par- 
ticle constant (no small mass particles are produced) and 
keeping the fractional mass of stars trapped in gas particle 
very small. We simply use an SPH-like gathering proce- 
dure. When the fraction of star-mass in a gas particle is 
higher than 5 %, we examine the star-mass content of its 
SPH neighbors. If there is enough star-mass in the neigh- 
boring particles to make one full star particle, the mass 



is gathered in the initial particle, which is then converted 
into a star. Using 15 SPH neighbors, the average star con- 
tent of a gas particle is smaller than 10%. This procedure 
is also applied to convert warm gas into cold gas, and stars 
into warm gas. Consequently we have to keep track of a 
"smoothing length" for all three types of particles. 

4. Formation of a stellar disk: reference case 

4-.1- Initial conditions 

The dark matter is distributed in a halo with a Plummer 
density profile: 



Pdm{r) = - 3-(l + ^) ^' 



We set Tdm = 5kpc, and Md,„ = 1.71O"M0. We in- 
troduce a cut-off at brdm- 

The baryonic matter (warm gas at 1=0) is dis- 
tributed in a rotationally supported thin disk defined by 
a Mi yamoto-Nagai density profile (Miyamoto & Nagai 
1975|) : 



+ {Tbm + z)2]5/2z3 



where, 



bm ' 



riim is the typical radius of the disk and Zdm is the typical 
vertical thickness. We use rf,m = 5 kpc and Zbm = 0.3 kpc, 
with cutoffs at 3rbm and 3zbm- The total baryonic mass is 
Mbm = 5.7 10^°Mq (1/3 of the dark matter halo mass). 
The angular speed ^(r) for the baryonic particle is derived 
from the density profiles of dark matter and baryonic mat- 
ter for global rotational support. Velocity dispersions are 
derived from disk stability criterion: 



3.36GE 



2n 



Vc^b 



where S is the surface density of the disk, and k the 
epicyclic frequency. Exact values are chosen to produce 
a Toomre stability criterion Q ^ 5 in the center, falling 
off to Q ~ 2 at the edge of the disk. 

Initially, all baryonic particles are SPH gas particle 
with a temperature of 15000K, and metallicity Z = 0. 

For the reference simulation we use the following pa- 
rameters. A Schmidt law with exponent 1.5 is used for star 
formation. The number of baryonic particle is 50000, the 
number of dark matter particle is 10000. This makes dark 
matter particles 15 times heavier than baryonic particles. 
The gravitational smoothing is 30 pc for baryonic matter 
and 300 pc for dark matter. The minimal SPH smoothing 
is accordingly set to 30 pc. The time step is 1 Myr and 
the integration is carried out over 2 Gyr. 
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Fig. 3. Reference simulation. Configuration of the three baryonic phases after 200 Myr (top) and after 1.2 Gyr 
(bottom) . 



4-- 2. results 

Face-on views of the various baryonic phases are presented 
in Fig. 3. The snapshots are taken at 0.2 Gyr and 1.2 Gyr. 
They show the formation of the stellar disk and the deple- 
tion of the gas phases. One interesting point in the stellar 
disk formation is that stellar clumps form in the outer, 
most unstable regions of the disk during the first few hun- 
dred Myr. They subsequently merge in the central part of 
the disk. We have encountered this phenomenon for a wide 
choice of parameters. Only with very high initial values of 
Toomre criterion {Q ^ 8) throughout the disk does it dis- 
appear. It is physically unlikely that such a Q = 8 disk is 
formed before any star formation takes place, since dissi- 
pation tends to decr ease Toomre criterion in the gaseous 
disk. Noguchi (199£) also forms clumps for a cold coUi- 
sional disk of gas and suggests that it can be a mechanism 
for the formation of a bulge. 

Another morphological aspect worth mentioning is the 
behavior of the spiral structure. It is strong in the three 
phases during the first few 100 Myr. It then fades, espe- 
cially in the stellar disk. Possible causes for this fading are 
the absence of external perturbations to drive the spiral 
structure, and the lack of resolution leading to unphysical 
heating of various baryonic components (see also Sec. 5) 



Fig. 4 shows the time-evolution of the mass of each 
baryonic phase. A large part of the warm gas is almost in- 
stantaneously turned into cold g thermal equilibrium 
is reached within a few time steps. Then mass exchanges 
proceed smoothly. Within 500 Myr, the warm gas mass 
reaches a quasi-constant value, as the balance between 
the gain from stellar mass loss and the loss from cooling 
into cold gas is reached. This is a self-regulated equilib- 
rium: if the warm gas mass is depleted through cooling, 
its density falls and the cooling becomes less efficient. If 
it increases through stellar mass loss, the density rises, 
the cooling becomes more efficient, increasing the rate of 
transfer to the cold gas phase. On long time scales (sev- 
eral Gyr), the warm gas mass would eventually decrease 
as stars grow older and release less warm gas. The cold 
gas phase is continuously depleted although the process is 
slow after 2 Gyr. We will examine in Sec 6 the influence of 
various parameters on the mass evolution of the phases. 
Fig. 5 shows the evolution of the star formation rate. The 
rate is high during the first Gyr, during which 75% of the 
gas is turned into stars. It decreases steadily with time. 
This is not fully into accord with the observations, which 
suggest more constant star formation rate. The missing 
factor here is probably the accretion of gas, which would 
replenish the gas disk and sustain star formation rate. 
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Fig. 4. Evolution of the composition of the baryonic mat- 
ter. The mass fraction of each of the three components is 
plotted against time. 

The surface density profiles of each baryonic phase at 
2 Gyr are plotted on Fig. 6. The stellar disk exhibits an 
exponential profile with a 2 kpc typical radius. The warm 
gas disk profile is also reasonably fitted by an exponential, 
with a rather larger typical radius of ~ 10 kpc. The cold 
gas disk falls in between. Finally, Fig. 7 shows that the 
circular velocity profile exhibit little evolution during the 
first 2 Gyr of the simulation. The profile is rather flat at 
r > 5 kpc, with a typical velocity of ~ 250 km s~^. 

4-3. Modification of the parameters 

The relative complexity of the model that we propose 
makes it impossible to study the action of all the param- 
eters independently. Some parameters have been tuned 
to produce physically acceptable results, and remain un- 
changed in our study. This is the case, for example, of the 
rate of dissipation through inelastic collisions in the cold 
gas: it is tuned so that cold gas looses a few percents of 
its velocity dispersion every 100 Myr. Another example 
is the temperature of the phase transition between warm 
and cold gas, set at 11000 K. Changing the value would tilt 
the mass-equilibrium between the phases. We have tuned 
it in connection with other parameters such as the heating 
rate F to produce similar masses for the two components 
after 2 Gyr in Model 1. 

The modifications we have chosen to study define 10 
models described in Fig. 8. Model 1 is the reference model 
described in Sec. 4. In model 2, heavy dark matter parti- 
cles have the same gravitational softening as light baryonic 
particles. In model 3, the dark matter halo is not realized 




0.5 1 1.5 2 



Time (in Gyr) 

Fig. 5. Star formation rate in the reference simulation, 
as a function of time. The star formation law is a Schmidt 
law with exponent 1.5. 

with particles, but using the static analytic profile given in 
Sec. 4.2. Model 4 includes 400 000 baryonic particles and 
100 000 dark matter particles with a unique 15 pc grav- 
itational softening, the time step is reduced to 0.5 Myr. 
In model 5, a Schmidt law with exponent 1 instead of 1.5 
is used to compute the star formation rate. In model 6, 
the stellar mass loss process is simply shut off, and an 
instantaneous recycling approximation is used for metal 
enrichment. In model 7, supernovae feedback is switched 
off. In model 8, the kinetic part of the supernovae feed- 
back is switched off. In model 9, the supernovae feedback 
is purely thermal again, but gas particle are heated up to 
100 000 K instead of 20 000 K. Last, in model 10, the mass 
of the dark matter halo is three times smaller than in the 
reference model. 

5. Disk thickness: the numerical heating 

Numerical simulations of galactic disks usually focus on 

the baryonic matter behavior which can be easily com- 
pared to observations. It is then tempting for a given CPU 
cost to increase the scale resolution of baryonic matter at 
the expenses of the dark matter, using many light baryonic 
matter particles and few heavy dark matter particles. This 
is the case in most of the simulations in this work, where 
the mass ratio between dark matter and baryonic matter 
particles is between 5 and 15. This raises the question of 
how to handle two-body relaxation within each compo- 
nent, and more importantly between the 2 components. 
The two-body relaxation time-scale is usually controlled 
by the value of the gravitational smoothing length. There 
is actually an optimal value which limits both two-body 
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Differences with model 1 


Model 1 


Reference model (see Sec. 4) 


Model 2 


Unique gravitational softening (30 pc) 


Model 3 


Analytical dark matter halo 


Model 4 


High resolution simulation: 500 000 particles 


Model 5 


Schmidt law with exponent 1. 


Model 6 


No stellar massloss 


Model 7 


No supernovae feedback 


Model 8 


Purely thermal supernovae feedback 


Model 9 


High thermal supernovae feedback 


Model 10 


Lighter dark matter halo 



Table 1. Description of the various models used in the simulations. 



relaxation and the bias in the force evaluation (see, e.g. 
Romeo 1997, Athanassoula 200C). How to deal with sev- 



eral particles species with different masses is not as well 
established. The main lead is to use variable softening pa- 
rameters. Dehnen ( |2001| ) presents an analytic study of the 
effect of using variable smoothing length on th e bias in 
the force computation. Binney & Knebe ( 2001 ) present 
a numerical study including two species of particles with 
different masses in a cosmological context, and find that 
using different softening values can reduce mass segrega- 
tion, an effect arising from two-body relaxation between 
the different species. 

Our early simulations showed a thickening of the stellar 
disk, ft was quickly realized that the heavy dark matter 
particles were responsible for the heating of the baryon 
disk. Being organized in a spherical halo, they have a 
large velocity dispersion perpendicular to the disk which 
is transmitted to the baryonic particles. We have stud- 
ied the importance of this effect under various simulation 
settings. 

The thickening of the stellar disk was compared for 
models 1, 2, 3 and 4. Fig. 9 shows the aspect of the disk 
seen edge on after 1.8 Gyr. Model 2, with a constant soft- 
ening length of 30 pc produces the thickest disk. Model 
1, with softening 30 pc for the baryonic matter and 300 
pc for the dark matter, produces a thiner disk, but still 
quite thicker than models 3 and 4. Model 3, with an ana- 
lytic halo component removes all heating of the disk due 
to the dark matter particles, thus providing a reference 
point. Interestingly model 4, including 500 000 particles 
also produces a very thin disk. The relaxation time has 



increased enough with the number of particles to produce 
little thickening within 2 Gyr. 

Fig. 10 gives a quantitative evaluation of the phe- 
nomenon. We plot the average vertical velocity disper- 
sion of the baryonic matter within a radius of 7 kpc as a 
function of time for the four models described above. Ini- 
tially the matter is mainly in the form of dissipative gas, 
which explains the initial decrease of the velocity disper- 
sion. Conclusions pertaining to the heating by dark matter 
can be drawn from the latter times ( > 0.5 Gyr), when the 
disk has settled into a quiet evolution. The conclusions are 
the same as drawn from Fig. 9. The only new element is 
that, in model 4, the heating is not completely suppressed, 
which was not obvious in Fig. 9. 

Our main conclusion is that, while different smoothing 
length helps reduce the transfer of vertical momentum be- 
tween dark matter and baryonic matter, a greater number 
of particle is the best solution. 



6. Sensitivity of physical quantities to parameter 
variations 

6.1. Mass equilibrium between the phases 

As explained in Sec. 4.4 and shown on Fig. 4, the param- 
eters of the various physical processes have been chosen 
in the reference simulation to produce a mass equilibrium 
between the phases that is only slowly evolving after 2 
Gyr. The gas fraction of the baryonic mass is then about 
21 % and slowly decreasing, which is compatible with the 
value of ^ 10% observed in present days, ^ 10 Gyr old 
galaxies. 
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Fig. 6. Surface density profiles of the various baryonic 
components of the disc as f mictions of the radius. The 
density scale is logarithmic. The profile are roughly expo- 
nential with 10 kpc typical radius for the warm gas, and 
2 kpc for the stars. 



Most of the physical processes included in the model 
can affect the mass balance between the phases. We have 
studied the effect of varying some of the parameters in- 
volved in the processes modeling. The result are summa- 
rized in Fig. 11. The mass fraction of each baryonic com- 
ponent is given at 0.5 Gyr , which marks the end of the 
initial most unstable period., and at 2. Gyr, the end of the 
simulations. 

Model 5 shows the effect of changing the exponent from 
1.5 to 1. in the Schmidt law. The constant C in eq. 3 is 
adjusted to yield similar overall star production; Model 5 
produces less stars in the initial period of high gas den- 
sity, and more in the later period of low gas density. The 
final fraction of cold gas is smaller in Model 5 than in the 
reference model, more of it having been turned into stars. 
The warm gas fraction is barely affected. 

Model 6 is a run where stellar mass loss has been 
switched off. It is interesting to notice that, although stars 
eject warm gas particle, it is the cold phase that is mostly 
depleted when mass loss is switched off. We interpret this 
as follows. The warm gas equilibrium mass which appears 
on Fig. 4, is associated with a critical density fixed by the 
cooling process. If the gas falls below this density, trans- 
fer to the cold gas phase almost stops. The injection of 
mass from stellar mass loss docs barely increases the den- 
sity above the critical density since the cooling increases 
in proportion to the density. Its main effect is actually to 
maintain a mass-flow from stars to warm gas to cold gas. 
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Fig. 7. Rotation curve of the disc (circular velocity as a 
function of radius) at the beginning and the end of the 
reference simulation. The curve shows little evolution and 
is rather flat in the outer part of the disk. The dark matter 
contribution is large, producing a stable disk. 



If this flow dies up at the source, by switching off stel- 
lar massloss, the warm gas settles at its critical density 
and stops producing cold gas. Cold gas is then steadily 
depleted through star formation. 

Models 7, 8 and 9 study the effect of supernovae feed- 
back. From the result of model 7, where SN feedback is 
absent, we assert that SN feedback does not affect the av- 
erage star production, but that it has a strong influence 
on the balance between cold and warm gas. We can guess 
that the evaporation of cold gas by thermal input from 
SN, associated with the damping of thermal evolution is 
a more important factor in this regard than the kinetic 
feedback. This is checked in model 8, where SN feedback 
is purely thermal. The fraction of each baryonic compo- 
nent is very similar to the reference model. Model 9 shows 
that using a high thermal feedback does not change the 
balance much more. 

The largest variation in this study is from 5.1 % to 
17.6 % of warm gas at 2 Gyr between models 6 and 7. 
The largest variation with respect to the reference model 
is from 17.6 % to 10.2 % of warm gas at 2 Gyr. These 
are limited variations which point to the stability of the 
model as far as the composition of the baryonic matter is 
concerned. 



6.2. Metallicity profiles 

We establish metallicity profiles at 2 Gyr to test the reli- 
ability of the overall star formation process. Let us recall 
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Model 


At 0.5 Gyr 


At 2. Gyr 




Stars 


Warm gas 


Cold Gas 


Stars 


Warm gas 


Cold gas 


1 


58.7 % 


12.5 % 


28.8 % 


78.9 % 


10.9 % 


10.2 % 


5 


53.7 % 


13.3 % 


33% 


84.8 % 


10.3 % 


4.9% 


6 


65.4 % 


11.8 % 


22.8 % 


86.2 % 


8.7% 


5.1 % 


7 


58% 


7.6% 


34.4 % 


77.5 % 


4.9 % 


17.6 % 


8 


61.2 % 


11.3 % 


27.5 % 


79.5 % 


10.4 % 


10.1 % 


9 


61.4 % 


13.2 % 


25.4 % 


78.5 % 


11.2 % 


10.3 % 



Table 2. Composition of the baryonic matter at 500 Myr and 2 Gyr in different models. See Fig. 8 for properties of 
the models. 




Model 3 Model 4 



Fig. 8. Edge-on views of the stellar component of the disk 
of different models, 1.8 Gyr after identical initial configu- 
rations. Models exhibit different degree of thickening. 

that metal enrichment is included in the simulation only 
as a diagnosis, without any influence on the cooling of the 
gas (see Sec. 2.2.1 for the justification). The metallicity 
profiles of the gas component for models 1, 5 and 6 are 
presented in Fig. 12. All three roughly fit an exponen- 
tial profile. As can be expected the metallicity gradient is 
steeper in models 1 and 6, with a Schmidt law of expo- 



35 I ' ' ' ' ' ' ^ 

0.5 1 1.5 

Time (in Gyr) 

Fig. 9. The vertical velocity dispersion of the stellar disc 

is averaged in the inner part of the disc ( r < 7 kpc). 
Its evolution is plotted for different models, tracking the 
thickening of the disk. 



nent 1.5, than in model 5, with a Schmidt law of exponent 
1. In Model 6, stellar mass loss is switched off, and met- 
als arc rctiirncd to the gas using the instantaneous recy- 
cling approximation. For every unit mass of star formed, 
0.4 unit mass of gas is instantaneously enriched using a 
yield y = 0.02. The difference with the reference simula- 
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Fig. 10. Metallicity (log scale) as a function of radius 
for different models (see Fig. 8). All profiles are roughly 
exponential. 



tion shows at small radii, where instantaneous recycling 
produces a steeper metallicity gradient. Indeed this is the 
region where star formation is the most active, and the in- 
stantaneous recycling, which does not allow enriched gas 
to be carried out of star forming regions, may overestimate 
the process of enrichment. 



6.3. Bar formation in a light dark matter halo 

In the reference simulation, we use a dark matter halo 
three times heavier than the baryonic disk within a 15 
kpc radius. This kind of configuration produces a stable 
disk and a rather quiet evolution which are convenient to 
test the model. However, it prevents the formation in the 
disk of bars and grand design spiral structures, in partic- 
ular m = 2 modes. Such features appear when the disk is 
less stable, that is for example if the dark matter halo is 
lighter or less concentrated. Moreover, disk rotation curves 
derived from observations tend to show that the contribu- 
tion of the dark matter to the circular velocity in the inner 
part of the optical d isk is small for typical spiral galaxies 
(e.g. Sofue & Rubin |200l| ). 

In model 10, we set the mass of the dark matter halo 
equal to the mass of the baryonic disk, that is three times 
lighter than in the reference simulation. All other param- 
eters are unchanged. The circular velocity of the disk in 
the outer flat part of the rotation curve is ^ 200 km s^^ 
in model 10, and ~ 250 km s~^ in the reference simula- 
tion. Fig. 13 shows the relative contribution of the dark 
matter halo to the circular velocity in the disk at t = 0, 
for models 1 and 10. The quantity plotted is: 
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Fig. 11. Relative contribution of the dark matter halo to 
the disc circular velocity as a function of radius (see main 
body for exact definition). In model 10 (light dark matter 
halo) , the inner part of the disc is dominated by baryonic 
matter. 



v(Pdm) 

v(Aot) 

the ratio between the circular velocity which would result 
from dark matter alone and the actual circular velocity. 
We can check that, in the inner part of the disk, model 10 
is closer to reality, with a small relative contribution from 
the dark matter. 

The analysis of the evolution of the baryonic matter 
composition doesn't show any qualitative difference be- 
tween model 1 and 10. Fig 14 shows the configuration of 
the various components at i = 900 Myr, for the reference 
model and for model 10. The existence of a bar and a 
m — 2 spiral mode is obvious in model 10, especially for 
the warm gas component. Simulation plots show that the 
bar is actually present during several 100 Myr. 

7. Conclusions 

The first motivation for this work is the recognition that 
the cold gas present in the ISM has morphological and dy- 
namical properties very different from those of the warm 
diffuse gas. Semelin and Combes (2000) and Semehn and 
Combes (2002) explore some of the properties of the cold 
gas and its effect on the galaxy formation process. Those 
were local studies, where the small scale fractal structure 
of the gas is fully taken into account. In the present work 
we try, using a simple model, to include cold gas physics 
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Fig. 12. Face-on plots of the four type of matter at 900 Myr for the reference model (top) and for model 10 (bottom). 



in a global model of galaxy formation and evolution. Ide- 
ally, what we learn from local studies of cold gas physics 
could be included as " sub- resolution" processes in global 
simulations. This interplay will hopefully develop in the 
future. 

This work presents a multiphase numerical model for 
galaxy formation and evolution. The model takes into ac- 
count four different types of matter obeying difTcrcnt dy- 
namics: dark matter and stars, warm gas and cold gas. 
Dark matter and star particles obey gravity only, com- 
puted in a tree algorithm. Warm gas particles follow self- 
gravitating hydrodynamics implemented by a tree-SPH al- 
gorithm, and cold gas obeys gravity and undergo inelastic 
collisions (sticky particle scheme). Warm and cold gas are 
treated as two different fluids, and this is one of the inno- 
vations in the present work. In galaxies, gas and stars ex- 
change matter and energy through a number of processes. 
The model takes into account the following: heating and 
cooling of the gas, star formation, stellar mass loss, ki- 
netic and thermal feedback from SN explosions, metallic 
enrichment. We use a continuotis massloss model which 
goes beyond the usual instantaneous recycling approxi- 
mation and has not been integrated in such a multiphase 
model before. 

Sec. 4 follows the evolution over 2 Gyr of an initially 
pure gas disk within a dark matter halo. The resulting 
galaxy exhibits properties similar to those of observed 
Milky Way like galaxies. We present in particular the evo- 
lution of the composition of baryonic matter, the star 
formation history, the surface density profile of the var- 
ious components and the initial and final rotation curves. 



These quantities (except for star formation history) is rea- 
sonably similar to observed values. 

The formation of a stellar disk from a pure gas compo- 
nent gives rise to a large number of self-gravitating clumps; 
those heat the disk, merge, and are driven towards the cen- 
ter through dynamical friction, where they contribute to 
form a bulge. The cold gas phase is always more concen- 
trated than the warm gas phase. Although the continuous 
stellar mass loss replenishes the disk in gas along its evolu- 
tion, the star formation rate decreases exponentially with 
a short characteristic time-scale of the order of 1 Gyr, 
and therefore gas infall would be required to reproduce 
the observed star formation history in giant galaxies like 
the Milky Way. 

Sec. 5 details the problem of the thickening of the stel- 
lar disk during the 2 Gyr of the simulation. This thickening 
is due to the two-body relaxation between star and dark 
matter particles. It is shown that using different, larger, 
smoothing length for the heavier dark matter particles 
reduces the transfer of velocity dispersion, but a better 
result is obtained by using an analytic halo or by increas- 
ing the number of particles to 500 000. We emphasize this 
point since many simulations of galaxy evolution imple- 
menting complex models (including ours) are still limited 
to a few 10^ particles for systematic studies. 

Sec. 6 explores the effect of several modifications to 
the model. We first study how the evolution of the bary- 
onic matter composition is modified under several type 
of circumstances, such as a different star formation law, 
different stellar massloss, or different SN feedback. The 
model is rather stable in this respect. The mass compo- 
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sition responds to these changes, but within reasonable 
limits. Constraining the final baryonic matter composi- 
tion from observations allows to discriminate between the 
various possible processes. The metallicity profiles are also 
computed for several models and are found to be exponen- 
tial, which validates the star formation process. 

If we try to estimate which of the processes and pa- 
rameters studied are the most significant for the overall 
evolution of the galaxy, it appears that taking into consid- 
eration a proper massloss scheme for the stars comes first. 
Indeed it affects noticeably the composition of baryonic 
matter, it affects the metallicity profiles, and, although it 
was not studied in detail in this work, it may affect the 
morphology of the galaxy by changing the sites of star 
formation. Next, comes the specific star formation law, 
which is of course central to reproduce the star formation 
history. And third, the existence of SN feedback which al- 
low to sustain a warm gas phase. It is an interesting point 
that the specific strength of the feedback is not a sensitive 
parameter, due to the cooling properties of the gas. Accre- 
tion of intergalactic matter, which was not studied in this 
work, may prove a key ingredient to reproduce the star 
formation history of spiral galaxies. It will be included in 
a future work. 

Finally a simulation was run with a lighter dark matter 
halo, that is a less stable initial configuration of the gas 
disk. As expected, the disc evolved to form a bar and an 
m — 2 spiral structure. 

Our model takes a large range of physical processes 
into account. For some of them, such as heating of the 
gas by the UV background, we use a very simple algo- 
rithm. There is room for improvement at this level. But 
more importantly, some of the physical processes, the star 
formation law in particular, are still not well understood. 
More simulation work is needed to validate or reject the 
various criterion in use. Finally, there is the important is- 
sue of modeling the cold gas dynamics. We used a simple 
sticky p article scheme. Other choices (such as Andersen & 
Burkert |2000|) are worth investigating in association with 
the SPH warm gas dynamics. Another possibility is to use 
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